interneurons <- readRDS("../../data/interneurons_kb.rds")
DefaultAssay(interneurons) <- "RNA"

interneurons = FindVariableFeatures(interneurons, selection.method = "vst", nfeatures = 10000)
all.genes = rownames(interneurons)
interneurons = ScaleData(interneurons)

interneurons = RunPCA(interneurons, features = VariableFeatures(object = interneurons))


interneurons = RunUMAP(interneurons, reduction = "pca", dims = 1:25)
interneurons= FindNeighbors(interneurons, reduction = "pca", dims = 1:25)


interneurons = FindClusters(interneurons, resolution = .9)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 19414
## Number of edges: 639483
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8362
## Number of communities: 23
## Elapsed time: 2 seconds
DimPlot(interneurons, reduction = "umap", label = TRUE, repel = TRUE,pt.size = 1)

plot_colors = read_delim("../../data/palettes/interneuron_timepoint_pal.txt",
                         delim = "/n",col_names = F) %>% pull(X1)

time_points <- unique(interneurons@meta.data["timepoint"]) %>% pull(timepoint) 

col_n <-set_names(plot_colors,levels(time_points))

#,order = rev(levels(time_points))
DimPlot(interneurons,
        reduction = "umap",
        repel = TRUE,
        group.by = "timepoint",
        pt.size = 1,
        cols = col_n,order = rev(levels(time_points)))+
  ggtitle(label = "")+ theme(legend.position = "none")

pre.critical <- levels(time_points)[1:5]

# Add new timepoint column for pre and post critical period
interneurons@meta.data <- interneurons@meta.data %>%
  mutate(timepoint2 = ifelse(timepoint %in% pre.critical,
                             "E14-P5",
                             "P7-Adult"))


DimPlot(interneurons, reduction = "umap", label = TRUE, repel = TRUE,pt.size = 1, group.by =  "timepoint2")

DimPlot(interneurons, reduction = "umap", label = TRUE, repel = TRUE, split.by =  "timepoint2",pt.size = 1)

saveRDS(interneurons,"../../data/interneurons_labeled.RDS")
DefaultAssay(interneurons) <- "RNA"

all.markers <- FindAllMarkers(interneurons,only.pos = TRUE,logfc.threshold = 0.25)

sig.markers <- all.markers %>%
  filter(p_val_adj < 0.05) 


sig.markers %>%
  filter(cluster == 6) %>%
  arrange(desc(avg_log2FC)) %>%
  head(20) %>%
  kbl() %>%
  kable_minimal()
p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
Adarb2 0 13.983543 0.824 0.530 0 6 Adarb2
mt-Rnr2 0 11.129866 0.970 0.854 0 6 mt-Rnr2
Hnrnpa2b1 0 8.444997 0.942 0.807 0 6 Hnrnpa2b1
Rpl21 0 7.708343 0.856 0.392 0 6 Rpl21
Tpt1 0 5.285972 0.943 0.824 0 6 Tpt1
Filip1l 0 4.589018 0.653 0.127 0 6 Filip1l
Rpl81 0 4.502508 0.957 0.842 0 6 Rpl8
Cpne4 0 4.200330 0.289 0.154 0 6 Cpne4
Atrx 0 3.748112 0.868 0.594 0 6 Atrx
Ccdc88a 0 3.693222 0.886 0.649 0 6 Ccdc88a
Fubp1 0 3.588860 0.797 0.655 0 6 Fubp1
Rps241 0 3.502384 0.965 0.854 0 6 Rps24
Srrm2 0 3.416876 0.936 0.716 0 6 Srrm2
Rpl39 0 3.393382 0.905 0.798 0 6 Rpl39
Thra 0 3.275159 0.780 0.543 0 6 Thra
Rpsa-ps10 0 3.041653 0.371 0.083 0 6 Rpsa-ps10
Dync1i11 0 3.034032 0.603 0.230 0 6 Dync1i1
Dclk11 0 2.961476 0.669 0.412 0 6 Dclk1
Sfpq1 0 2.950321 0.847 0.589 0 6 Sfpq
Nedd4 0 2.869791 0.824 0.476 0 6 Nedd4
write_csv(sig.markers, file = "sig_markers_1.csv")
# genes in the P7 clusters that have a log2 fold change >= 1
P7_genes  = sig.markers %>%
  filter(cluster ==6 & avg_log2FC>=1) 

# Convert gene symbols to ENSEMBL
geneconvert = bitr(rownames(P7_genes), fromType="SYMBOL", toType="ENSEMBL", OrgDb="org.Mm.eg.db")

# Convert to MGI IDs
mgi_ids <- select(org.Mm.eg.db, geneconvert$ENSEMBL, "MGI","ENSEMBL") %>% pull(MGI)

# Write MGI IDs for genewalk analysis
mgi_ids %>% as_tibble() %>% write_delim(file = "P7_genes.txt",col_names = FALSE)

go_enrich <- enrichGO(gene = geneconvert$ENSEMBL,
                      OrgDb = org.Mm.eg.db,
                      keyType = "ENSEMBL",
                      pvalueCutoff = 0.05,
                      ont = "BP",
                      pAdjustMethod ="BH",readable = TRUE)



# clusterProfiler::cnetplot(go_enrich,
#                           showCategory = 15,
#                           layout = "gem",
#                           cex_gene=2,
#                           cex_category=2)

go_plot <- barplot(go_enrich,showCategory = 15,font.size = 24)
  #scale_y_discrete(expand=c(0.2, 1))
go_plot$data <- go_plot$data %>% tail(-2)

go_plot

# dotplot(go_enrich,
#         split="ONTOLOGY",showCategory=15)+
#   facet_grid(ONTOLOGY ~ .,scales = "free")+
# scale_y_discrete(expand=c(0.1, 0))
table(Idents(interneurons)) %>% kbl(col.names = c("Cluster","Cell Count")) %>% kable_minimal()
Cluster Cell Count
0 2104
1 1882
2 1754
3 1442
4 1348
5 1285
6 1206
7 1149
8 1009
9 894
10 857
11 838
12 723
13 659
14 572
15 492
16 401
17 273
18 256
19 112
20 93
21 35
22 30
DimHeatmap(interneurons, dims = 1:15, cells = 500, balanced = TRUE)

DimPlot(interneurons, reduction = "pca", group.by = "timepoint",pt.size = 1)

DimPlot(interneurons, reduction = "pca", group.by = "timepoint",pt.size = 1,dims=c(2,3))

DimPlot(interneurons, reduction = "pca", group.by = "timepoint",pt.size = 1,dims=c(3,4))

DimPlot(interneurons, reduction = "pca", group.by = "timepoint",pt.size = 1,dims=c(4,5))

DimPlot(interneurons, reduction = "pca", group.by = "timepoint",pt.size = 1,dims=c(5,6))

DimPlot(interneurons, reduction = "pca", group.by = "timepoint",pt.size = 1,dims=c(5,7))

DimPlot(interneurons, reduction = "pca", group.by = "timepoint",pt.size = 1,dims=c(7,8))

interneurons[["percent.olf"]] <- PercentageFeatureSet(interneurons, pattern = "^Olfr")
interneurons[["percent.vmn"]] <- PercentageFeatureSet(interneurons, pattern = "^Vmn")

FeaturePlot(interneurons,features = "percent.olf",label = TRUE, repel = TRUE, min.cutoff = "q10", max.cutoff = "q90",pt.size = 1)

FeaturePlot(interneurons,features = "percent.vmn",label = TRUE, repel = TRUE, min.cutoff = "q10", max.cutoff = "q90",pt.size = 1)